home *** CD-ROM | disk | FTP | other *** search
/ Gold Medal Software 1 / Gold Medal Software Volume 1 (Gold Medal) (1994).iso / graphics / tierra40.arj / TIERRA / TRAND.C < prev    next >
C/C++ Source or Header  |  1992-09-09  |  3KB  |  102 lines
  1. /* trand.c  9-9-92  adapted from ``Numerical Recipes in C'' by Press,
  2. Flannery, Teukolsky and Vetterling.  Tdrand() returns a uniform random deviate
  3. between 0.0 and 1.0.  Tlrand() returns a uniform random deviate between 0 and
  4. 2^31 */
  5.  
  6. #ifndef lint
  7. static char     sccsid[] = "@(#)trand.c    1.5     7/21/92";
  8. #endif
  9.  
  10. #include "tierra.h"
  11. #include "extern.h"
  12.  
  13. #ifdef MEM_CHK
  14. #include <memcheck.h>
  15. #endif
  16.  
  17. #define M1 259200L
  18. #define IA1 7141L
  19. #define IC1 54773L
  20. #define RM1 (1.0/M1)
  21. #define M2 134456L
  22. #define IA2 8121L
  23. #define IC2 28411L
  24. #define RM2 (1.0/M2)
  25. #define M3 243000L
  26. #define IA3 4561L
  27. #define IC3 51349L
  28.  
  29. void tsrand(seed)
  30. I32s  seed;
  31. {   I16s  j;
  32.  
  33.     RandIx1 = (IC1 + seed) % M1;       /* Seed the first routine */
  34.     RandIx1 = (IA1 * RandIx1 + IC1) % M1;
  35.     RandIx2 = RandIx1 % M2;            /* and use it to seed the second */
  36.     RandIx1 = (IA1 * RandIx1 + IC1) % M1;
  37.     RandIx3 = RandIx1 % M3;            /* and third routines. */
  38.     for (j = 1; j <= 97; j++)    /* Fill the table with sequential uniform */
  39.     {   RandIx1 = (IA1 * RandIx1 + IC1) % M1; /* deviates generated by the */
  40.         RandIx2 = (IA2 * RandIx2 + IC2) % M2; /* first two routines. */
  41.         TrandArray[j] = (RandIx1 + RandIx2 * RM2) * RM1;
  42.     }                    /* Low- & high-order pieces combined here */
  43. }
  44.  
  45. double tdrand()
  46. {   double temp;
  47.     I16s j;
  48.  
  49.     RandIx1 = (IA1 * RandIx1 + IC1) % M1;
  50.                              /* Except when initializing, this is where we */
  51.     RandIx2 = (IA2 * RandIx2 + IC2) % M2;
  52.                              /* start.  Generate the next number for each */
  53.     RandIx3 = (IA3 * RandIx3 + IC3) % M3;                  /* sequence. */
  54.     j = 1 + ((97 * RandIx3) / M3);
  55.                        /* Use the third sequence to get an integer between */
  56.     if (j > 97 || j < 1) /* 1 and 97 */
  57.     {   FEError(-1200,EXIT,NOWRITE,"Tierra tdrand() This cannot happen.");
  58.     fprintf(stderr,"Tierra tdrand() TIERRA INTERNAL ERROR: %s",__FILE__);
  59.     exit(-1200);
  60.     }
  61.     temp = TrandArray[j];               /* Return that table entry, */
  62.     TrandArray[j] = (RandIx1 + RandIx2 * RM2) * RM1;  /* and refill it. */
  63.     return temp;
  64. }
  65.  
  66. /* ======================================================================= */
  67. #ifdef FUTURE
  68. float gasdev()
  69. {
  70.     static I32s iset=0;
  71.     static float gset;
  72.     float fac,r,v1,v2;
  73.  
  74.     if  (iset == 0) {
  75.         do {
  76.             v1=2.0*tdrand()-1.0;
  77.             v2=2.0*tdrand()-1.0;
  78.             r=v1*v1+v2*v2;
  79.         } while (r >= 1.0);
  80.         fac=sqrt(-2.0*log(r)/r);
  81.         gset=v1*fac;
  82.         iset=1;
  83.         return v2*fac;
  84.     } else {
  85.         iset=0;
  86.         return gset;
  87.     }
  88. }
  89. #endif
  90. /* ======================================================================= */
  91. #undef M1
  92. #undef IA1
  93. #undef IC1
  94. #undef RM1
  95. #undef M2
  96. #undef IA2
  97. #undef IC2
  98. #undef RM2
  99. #undef M3
  100. #undef IA3
  101. #undef IC3
  102.